Technical Report # 5-32688 
Contract Number NAS8-36955 
Delivery Order No. 122 


Turbulence Modelling of 
Flow Fields in Thrust Chambers 
(5-32688) 


Final Technical Report for the Period 
June 10, 1991 through September 13, 1992 


February 1993 


Prepared by: 


C. P. Chen * 

Y. M. Kim 
H. M. Shang 

Authors listed alphabetically 


Prepared for: 


NASA Marshall Space Flight Center 
Attn: Klause Gross 
Propulsion Laboratory 
Marshall Space Flight Center, AL 32512 


NASA 

National Aeronautical and 


Report Document Page 


|1 . Report No. 

Final Technical 


2. Government Accession No. 


ti. Recipient's Catalog No. 


W. Titile and Subtitle 


Turbulence Modelling of Flow Fields in Thrust Chambers 


p. Report Due 
September 1992 


|6. Performing Organization Code 
Research Institute, UAH 


y. Author(s) 

C. P. Chen , * Y. M. Kim, and H. M. Shang 
* Authors listed alphabetically 


p. Performing Organization Report No. 


p. Performing Organization Name and Address 

UAH Research Institute 
Rl E-47 

Huntsville, AL 35899 


|10. Work Unit No. 
5-32688 


|1 1. Contract or Grant No. 
NAS8-36955, D.O. 122 


[12. Sponsoring Agency Name and Address 

National Aeronautics and Space Administration 
Washington, D.C. 20546-001 
Marshall Space Flight Center, AL 3581 2 


[1 3. Type of report and Period covered 
Final Technical Report 

June 10, 1992 - Sept. 13, 1992 


fl 4. Sponsoring Agency Code 


p 5. Supplementary Notes 


1 6. Abstract 


Following the consensus of a workshop in Turbulence Modelling for Liquid Rocket Thrust Chambers, the current 
effort was undertaken to study the effects of second-order closure on the predictions of thermochemical flow 
fields. To reduce the instability and computational intensity of the full second-order Reynolds Stress Model, an 
Algebraic Stress Model (ASM) coupled with a two-layer near wall treatment was developed. Various test 
problems, including the compressible boundary layer with adiabatic and cooled walls, recirculating flows, swirling 
flows and the entire SSME nozzle flow were studied to assess the performance of the current model. Detailed 
calculations for the SSME exit wall flow around the nozzle manifold were executed. As to the overall flow 
predictions, the ASM removes another assumption for appropriate comparison with experimental data to account 
for the non-isotropic turbulence effects. 


17. Key Words (Suggested by Author(s)) 

18. Distribution Statement 

Turbulence, Combustion, Swirling Flows, 

TBD 

Non-Isotropic Effects, Compressible Effects 



9. Security Class, (of this report) [20. Security Class, (of this page) |21. No. of pages [22. Price 


Unclassified Unclassified I 24 + Appendices 


NASA FORM 1 626 OCT 86 






TABLE OF CONTENTS 


Abstract j 

Aeknowlegement 2 

1 . Introduction 3 

2. Physical Models 4 

1 . Eddy Viscosity Models 

2 . Second Order Models 7 

3. Numerical Implementation 13 

4. Results and Discussions 

5 . Conclusions and Recommendations \ 9 

References 20 

Appendix 1 . Publication Resulting From This Contract 23 

Appendix 2. Sample Inputs 24 

Figures 26 



ABSTRACT 


Following the consensus of a workshop in Turbulence Modeling for Liquid Rocket 
Thrust Chambers, the current effort was undertaken to study the effects of second-order 
closure on the predictions of thermochemical flow fields. To reduce the instability and 
computational intensity of the full second-order Reynolds Stress Model, an Algebraic 
Stress Model (ASM) coupled with a two-layer near wall treatment was developed. Various 
test problems, including the compressible boundary layer with adiabatic and cooled walls, 
recirculating flows, swirling flows and the entire SSME nozzle flow were studied to assess 
the performance of the current model. Detailed calculations for the SSME exit wall flow 
around the nozzle manifold were executed. As to the overall flow predictions, the ASM 
removes another assumption for appropriate comparison with experimental data, to account 
for the non-isotropic turbulence effects. 
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1. INTRODUCTION 

Thermochemical flow fields of liquid rocket thrust chambers are highly irregular in 
nature. The momentum fluxes as well as the scalar fluxes, due to the random velocity 
fluctuations, are usually much greater than the molecular (or laminar) fluxes. Proper 
descriptions of the turbulent mixing is the key of understanding and predicting reacting 
flow fields in thrust chambers. 

One of the most important features of the flow field analysis in the thrust chamber is to 
describe the turbulence structure, which correctly includes the variable density effects. In 
combusting flows, they are caused mainly by chemical reaction, mixture of gases with 
different densities, including multi-phase interactions, and strong distortion caused by 
shock boundary layer interactions. Complexity due to compressible combusting turbulence 
concerning the fluctuating density, leads to more non-linear equations to be solved and to 
more intractable modeling problems. 

Currently, methods of simulating turbulent flows can be classified in the following 
way. 1) Empirical Correlations, 2) Integral Methods, 3) One-point Closure Methods, 4) 
Two-point Closure Methods, 5 ) Large eddy simulation (L.E.S) and 6 ) Direct numerical 
simulation (D.N.S). Moving downwards in the list, each method requires more 
computational resources and fewer modeling assumptions. However, the higher level 
simulations also require more complex input data to calibrate the models. Such information 
is very difficult to obtain and often not as accurately known. Therefore, a simulation at a 
given level is not always more accurate than simulations at lower levels, which applies 
especially to level 2 , 3 and perhaps 4 . Due to the tremendous cost of simulating even the 
simple flows with high level L.E.S and full direct simulations on present day computers, 
and considering also the trade-off between accuracy and computation time, the one-point 
closure methods seem to offer the best compromise for engineering technology applications 
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at the present time. The one-point closure methodology includes: (i) Gradient Transport 
Models and (ii) Reynolds Stress Models. In (i) the construction of an eddy viscosity can be 
further classified into : (a) equilibrium models, which are essentially mixing length type 
models; and (b) non-equilibrium models, which include one and two equation models, 
Reynolds Stress and Algebraic Stress models. 

Combustion can affect turbulent transport through the production of density 
variations, buoyancy, dilatation due to heat release, influence on molecular transport, 
instabilities and so on. These effects are not well understood now. By and large, empirical 
closures and model equations are carried over from the Reynolds-averaging procedure for 
constant density non-reacting flows to the density-weighted (or Favre- averaged) form for 
turbulent reacting flows. In doing so, physical interpretations of the individual terms can be 
clearer than for the case with unweighted averaging. However, extra density fluctuating 
coupled terms, appearing in the transport equations, require extensive modeling efforts. 
New experimental data is needed to validate the modeling of these new terms. 

The main purpose of this proposed study is to incorporate a model which contains 
sufficient flow physics, encountered in the thrust chamber, and yet is computationally 
effective to include chemistry/droplet/multiphase effects into a CFD code, to systematically 
evaluate/assess the performance of the model. In the following, the proposed methodology 
and proposed tasks will be described. 

2. PHYSICAL MODELS 

In compressible flows, the statistical description of turbulent flow fields is derived 
by expressing the dependent variables as the sum of mean and fluctuating quantities and 
then ensemble-averaging the instantaneous transport equations. Due to the variations and 
fluctuations of the fluid density, which are more likely to occur in reacting flows with large 
temperature differences, it is advantageous [1,2,3] to use mass- weighted (or density- 
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weighted) averaging procedures tor describing the mean flow features. In chemically 
reacting flows the instantaneous velocities and temperatures of the fluid are averaged with 
density-weighting, and density and pressure are merely ensemble-averaged without 
weighting. Thus, for example, 


(1) 


P=P + P' 


and 


also 


Qi =puj / p 
puj =0 


( 2 ) 

(3) 

(4) 


Ui = -p U; / p * 0 


(5) 


Equation (5) is important, because most models of compressible turbulent flows in 
current use neglect such terms or their divergence, which may introduce problems in 
chemical reacting flows. 


The resulting statistical equations for the mass, momentum, energy fields as well as 


species in turbulent flows are 



T~(PUi) = 0 
OXj 


( 6 ) 


dpfij , 
at 


a _ 


a - 


— <pa jDj ) = -_p + __ (MS .. ) - — (p Uj up 


3x: 


(7) 
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where the mean strain is 


3xj 3xj 


3xj 3xj 


(O k 


(9) 


and 


S u 


^ + i^_2 8 JK 

3xj 3xj 3 i] 3x k 


Note that the density-weighted averaging produces forms in the mean momentum 
equations which are similar term by term to their incompressible counterparts. As 
mentioned in the introduction section of this report, the main task of the current one-point 
statistical turbulence modeling is to evaluate the Reynolds stresses, which appear as the last 
term in equation (7), as well as second order terms in the energy and species equations. 

1). Eddv Viscosity Models 

At the simplest level of modeling, an explicit algebraic constitutive relationship 
between the Reynolds stress and the mean strain rate of the flow has been postulated. This 
requirement utilizes the concept of an eddy viscosity as a ratio of proportionality. From a 
numerical point of view, this relationship is very convenient, since the eddy viscosity can 
be combined with the fluid kinematic viscosity to form an “effective “ viscosity. It thus 
allows simple modification of a CFD code, originally developed for laminar flows, to 
account for turbulence. 
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The simplest constitutive relationship for expressing Reynolds stress is the classic 
Boussinesq form, which is still the most popular one: 




At this level of the model, the major concern is the estimation of the eddy viscosity p, 

In general, the eddy viscosity is expressed in terms of a characteristic turbulence velocity 
scale u, a characteristic length scale 1 , and an empirically determined constant c^. 

Turbulence models, based on the Boussinesq expression, include zero-equation (mixing 
length) models [5,6,7], two-equation (with or without modifications of variable density 
effects) [8-13], and multiple-scale models associated with eddy viscosity estimations. This 
level of modeling is still the most widely used model in most engineering calculations. 
However, the major drawback of the Boussinesq expression is that the constitutive 
assumption implies an isotropic turbulence field. Past results show that for many complex 
flows, especially for shear driven flows such as the ones experienced in the coaxial injector 
jet flows in the combustor of a thrust chamber, this type of modeling usually fails to predict 
the turbulent flow field correctly. 


Second Order Models 


Recently, a workshop was held [16] in which the currendy used turbulence models 
were reviewed. Based on the consensus of the workshop, models based on second-order 
closures were recommended. During the following activity, addressing the second order 
modeling level, the constitutive assumption such as eq.(10) was totally abandoned and a 
transport equation for the Reynolds stress was derived. By manipulating the instantaneous 
Navier-Stokes equation, the field equations representing Reynolds stress can be obtained: 
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( 11 ) 


D — - 

^pUiUj=Pij + D ij + jt ij + C ij -e ij 

To date, equation (11) still represents the most complex turbulence model, used in 
computational fluid dynamic problems, in which anisotropy and extra complex strains can 
be automatically accounted for. The symbols and their physical mechanism as well as their 
mathematical forms are given in Table 1 . These terms contain a multitude of new turbulence 
moments involving pressure fluctuation, velocity fluctuations, velocity gradients, and third 
order correlations. Except for the production terms, the other terms require appropriate 
modeling to close the set of equations. Detailed experimental data are required to guide 
these modelings. However, at this stage, there are no reliable data available especially for 
reacting shear layers to give any guidance on the model development. Consequently, the 
validity of several terms in the model can only be inferred indirectly from the behavior of 
the mean motions and some budget profiles across some simpler flows. For three- 
dimensional flow calculations, six Reynolds-stress equations have to be solved which 
introduce a severe computer cost penalty. To alleviate this problem, but retain the 
advantages of this level of modeling, a simplified version is recommended. 


The Algebraic Stress model proposed by Rodi[4] is based on the assumption that the 
difference of advection (convection + diffusion) of Reynolds stress and its diffusion 
transport is proportional to the corresponding difference in turbulent kinetic energy, which 
implies 


Thus 


D 7^ r~. P u i u j .Dk ^ x 
-pu iUj -D ij = -^(_- Dk ) 


pUjUj 


k (Pjj + Kjj -£jj + C,j) 


Pl- *+“ Cl 


( 12 ) 
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Terms in Reynolds Stress Equation 


Pi] = Production Tensor = 




du i 


dxi 


4 - u'j u 


n du ; ] 
k 8x it J 


D{j = Diffusion Tensor 
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[pw'u'ju'/' + <\ku''p' + Sjku'jp' - (pSiktPj t- 
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du*? 

(JXi 


Cij =Compressibilitv Tensor = -fu"# 2 - + u"#2-l 

l i oxj j ax t J 

6,7 = Dissipation Tensor = §j£ + 


Terms in Turbulent Kinetic Energy Equation 


T’fc = Production Term = 

1 i J dxj 

D k = Diffusion Term = -gf- -{pufk + <7) + gf '-(/ji&Juf) 
Ck = Compressibility Term = p'S^- - i£- 

e = Energy Dissipation = /45q 


Table 1. Mathematical Forms of Reynolds Stress and 
Kinetic Energy Equation 


where P k ,C k and e are the contractions of their corresponding terms. With this algebraic 
form the cost effectiveness of just using a kinetic energy equation, rather than the six 
Reynolds-Stress field equations, is reduced. The kinetic energy equation is obtained by 
contracting the Reynolds stress equation and the modeled e equation as follows: 

Dk „ 

— -P k +D k +C k -e ( 13 ) 

The Algebraic Stress model has considerable appeal, as it offers the advantage of the 
Reynolds Stress model that is, all effects that enter the transport equations for pujuj' 

through the source terms, such as body forces (rotations, streamline curvature, buoyancy), 
anisotropic strain field, high order compressibility, and wall damping influence, can be 
incorporated. On the other hand, instead of solving six extra partial difference equations, 
only six algebraic equations must be solved. The resulting model is much more reasonable 
for application than the Reynolds Stress model and eliminates some uncertainties, 
associated with assigning proper boundary conditions for the Reynolds stresses. However, 

this model relies on many of the closure assumptions used for Reynolds Stress Models. 
Hence, only by properly validating the closure involving and k 1} terms, can we 

have confidence in the simpler models which are derived from it. 

For non-equilibrium models, at least one of the characteristic scales was estimated by 
a transport equation to account for the history effects. In most cases, the characteristic 
turbulence velocity scale is estimated from taking the square root of the turbulent kinetic 
energy which is governed by equation (13). The characteristic length scale is then 
determined from a combination of the turbulent kinetic energy and the other scalar 
turbulence quantity. The turbulence kinetic energy dissipation rate e and the mean square 
root of vorticity to seem to be most popular choices for length-scale determining transport 
equations from many other choices[9,12]. The equation for e can also be derived from the 
original instantaneous Navier-Stokes equation. However, the more important issue is the 
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interpretation of e . By definition, the destruction of turbulent energy by viscous action 
occurs at the finest scale of the fluctuation and is locally isotropic. Using this quantity to 
determine the characteristic length scale for momentum mixing, which reflects essentially 
the large-eddy motions, relies heavily on the assumptions of the energy “cascade” 
process[8]. That is, the turbulence energy dissipation rate is controlled by the rate, at which 
energy cascades from large to small scale eddies. Thus, the “modeled” transport equation 
for e depends heavily on the analogy of the corresponding turbulent kinetic energy 
equation. The form most frequently used is 

K (pe) = ^ ( £^ ) + C >Pf P >< ~C 2 Py+ Compressibility Effects ( 15 ) 

The resultant transport equations are also highly sensitive to the model constants Q 
and C 2 . For example, a change of Cj by one percent would decrease the predicted 
speading rate of a subsonic round jet by as much as 5%. The calculation of the length-scale 
governing equation itself is probably the weakest link in this level of modeling and may 
ultimately force an adoption of multipoint methods. Other two-equation models are the 
k-o) 2 model of Wilcox and Rubesin[9], or the k 1/2 -aj model of Coakley[17] and 
Bardina[l 1], where 1/to is related to the characteristic eddy life time k/o. Experience with 
the two-equation models indicate the inability of all models to adequately predict the extent 
of separation caused by shock/boundary layer interaction. Some improvement is obtained 
by introducing a compressibility correction in the original incompressible version of the 
k - e model, simply by changing the constant that multiplied the dilation terms[18] or by 
making the model constant Mach number dependable on some ad hoc assumptions! 13], 

The major task of the Algebraic Stress Model is to model the unknown terms on the 
right hand side of equation(12). Specifically the pressure strain tensor rc, r which controls 
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the redistribution of turbulent energy among the normal stresses through the interaction of 
pressure and the strain rate, and the viscous dissipation tensor e ir The pressure strain term 
is modeled through three mechanisms: resulting from purely turbulence interactions 

known as “retum-to -isotropy”[19], n i}2 involving interactions between the mean strain rate 
and turbulence, the so-called “rapid term”, and tt ijw relating the effects of solid boundaries 
on both Tty! and jtjj 2 . 


In the present study, the most frequently used linear model of Rotta [19] is adopted 
for the retum-to-isotropy: 

reijl = ~ C, t (pUiUj ~3 ijk) (15) 

More complicated non-linear models, such as the models of [2,20], have been proposed, 
however, these have shown no significant improvement over Rotta’s model. The rapid 
term is approximated by the isotropization production (IP) model, suggested by 
Launder[21], thus: 

rtij2=-C 2 (PSj-|6 ij P lt ) (16) 

in which 

P ‘4 P ° 

Finally, the dissipation tensor is modeled by 

c _ 2 C 

e ‘j - 3 (17) 

where e is the turbulent energy dissipation rate that can be obtained from equation! 1 1). 


With these models, the final formulation of eq.(12) results in: 


P U iUj 

k 



i-c 2 

(C, - l)e + P k 



(18) 
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It should be noted that the left hand side of the equation is a non-dimensional measure 
of anisotropy of turbulence. 

For the term tt 1jw , various researchers[20,22,23] have developed models, 

accounting the echo effect”, to damp the normal velocity fluctuation and redistribute its 
energy to the other two turbulence intensities. Available models for * ijw involve many 

complicated terms, and it is not clear, how these formulations can be directly applied to 
complex geometries. To avoid uncertainties due to the boundary effect, we have adopted a 
composite modeling approach, which implicidy accounts for the near wall effects. In this 
study, a two-layer approach is implemented. In the fully turbulent region away from wall, 
the form of the ASM as described above is adopted. In the near wall region, including 
overlap regions and viscous sublayer, the one equation k - 1 model is used, with a scalar 
eddy viscosity in the inner layer, while the influence of extra rates on Reynolds stress 
distributions are accommodated in the fully turbulent region. The present approach 

provides improved resolution over the RSM/wall function approach, commonly adopted in 
the literature [22]. 


The matching point for this two layer ASM model is chosen as y + =200. Within the 
matching point, the Reynolds stress tensor is calculated from 


in which 


-p Uj uj = c„k V(f*. + - f «S fjj|4 - f V* 


*n =Cin[l-exp(- ^ Ct 2 ^ )] 
An A 


(19) 


( 20 ) 
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The near wall shear stress damping provided by 1^ implicitly accounts for near wall 
pressure strain, though it is only appropriately correlated to the viscous damping effects. 


Although the nonisotropic stress model of Eq.(18) is used to evaluate turbulent fluxes 
in the momentum and kinetic energy equations, the gradient diffusion approximation is 
retained for the turbulent heat flux terms and turbulent mass terms in the energy and species 
equations. The effective transport coefficients are modeled through the introduction of the 
turbulent Prandtl number and turbulent Schmidt number. Thus 


V p, 3h 

pu.h = 

Pr t dx { 


( 21 ) 


and 


pUjC = 

Sc t 3xj 


( 22 ) 


These turbulent Prandtl and Schmidt numbers are set to a constant value of 0.9 as 
commonly used in the literature [34], 


3. NUMERICAL IMPLEMENTATIONS 

The governing equations(6-9) are mapped into a general body fitted coordinate 
system by the standard transformation formulae[24]. The dependent variable <t> in the 
general coordinate system can be expressed in a compact form as follows: 


dt 




where u, V are the contravariant variables that represent convective fluxes 


U = gllU + gl2 v 
V = g 2 iu + g 22 v 
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and r<j, and S^are the associated diffusivity and source terms for the variable ® 
(-u,v,h,k,e, etc.). The detailed expressions of the source terms and the cross derivative 
terms due to grid non-orthogonality are available in [25], These equations are then 
discretized using a finite difference method, based on the control-volume formulation on a 
non-staggered grid arrangement for all dependent variables. The velocity-pressure coupling 
was resolved earlier by the PISOC algorithm in a time-marching fashion. The salient 
features of the current MAST-2D include: strong conservation form with Cartesian 
components as dependent variables, colocated grid/variable arrangement, pressure- based 
PISO-C algorithm, high order Chakravarthy-Osher TVD scheme, and conjugate gradient 
(CGS) matrix solver. The current numerical method is capable of computing internal and 
external flows that are laminar, turbulent, separated or attached, incompressible or 
compressible and requires no smoothing or explicit under-relaxation other than implied by 
variations in the time step. Table 2 provides an updated list of features and capabilities of 
the MAST-2D code. Further details on the numerical method are provided in [25]. 


Equation(18), which represents a nonlinear coupled system, is implemented in the 
flow solver by lagging the turbulent kinetic energy production terms P by one time-iteration 
step. The linearized system is then solved at every grid point directly by Gaussian 
elimination. The final model has been validated for several flows including a compressible 
developing boundary layer flow, a separated flow, and a swirling flow. A complete thrust 

chamber flow field calculation has been executed also. The computed results will be 
discussed below. 
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Table 2. Updated features 
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4. RESULTS AND DISCUSSIONS 
Boundar y Laver Flows 

The model, developed in this study, is used first to calculate compressible flat plate 
boundary layers on adiabatic as well as cooled walls, and the results are compared with 
benchmark experimental data of [26, 27]. The calculations were carried out over the Mach 
number range, 0<Ma<5, for the adiabatic wall condition and over the temperature range, 
0.2<Tw/Taw<l, for the cool wall case. Here, Ma is the free stream Mach number , Tw is 
the prescribed wall temperature, and Taw is the adiabatic wall temperature. Calculations 
were done by either direct integration to the wall, using the damping functions described 
above or integration to the inertia sublayer, using the wall function. In Figure 1, the 
normalized skin frictions, predicted by the ASM and the two-equation k - e model are 
compared with the Van Driest curve and experimental data for a range of external Mach 
numbers. The current ASM model gives even better predictions, compared to the recent 
compressibility-corrected k-e model predictions [28], and is also in good agreement with 
the full Reynolds Stress Closure results [30]. The predictions for the cooled wall are 
shown in Figure 2. As pointed out by [29], if the cooled- wall flows are to be predicted 
correctly, turbulent heat fluxes and their near- wall behavior need to be realistically 
modeled. The current two-layer approach appears to model these aspects satisfactorily. In 
Figure 3, the near wall mean velocity profiles are plotted in terms of wall parameters, based 
on friction velocity, wall shear stresses and molecular viscosity. The calculated profiles 
based on the ASM/two-layer model are in agreement with data over a wide range of y+. 
The slope of the calculated profiles are also roughly parallel to that determined from 
measurements. The predictions of the near wall flow can be further examined from the 
profiles of the normalized kinetic energy and the normalized Reynolds stress in Figures 4a 
and 4b. These results are consistent with the asymptotic analytical results deduced from the 
direct numerical simulation calculations [29]. 
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2.-D Backward- Facin g Step Flow 

The second test case is the two-dimensional backward-facing step flow, involving 
massive flow separation. The experimental data of Driver and Seegmiller [3 1 ] was used for 
comparison. The treatment of this problem follows closely the procedure described in [25], 
Thus, only the relevant second order results are presented here. The predicted values of the 
reattachment length of the recirculation zone behind the step, using the k - e model and the 
ASM closure, are 4.76 and 5.94 step heights respectively. The reported experimental data 
on the reattachment length is about 6.1 with some fluctuations of about a 0.2 step height. 
The current ASM model performs much better than the eddy-viscosity k - e model for 
recirculating flows. In terms of turbulence quantities, comparisons of the predicted 
streamwise and lateral turbulence intensities, as well as the Reynolds stress component 
obtained from the solutions of ASM Eq. (18), and estimates from the constitutive equation 
(10) by using the k - e model, are shown in Figure 5. The results were normalized with 
the inlet centerline velocity. The non-isotropic imbalance of the turbulence intensities as 
well as the level of Reynolds stress are much better predicted by the ASM closure. 

Confined Swirling Jet Flow 

The next case involves flows with separations and streamline curvature due to swirl. The 

confined swirling jet, experimentally carried out by Roebuck and Johnson [32], is chosen 

for the model assessment study. The swirl number, S, obtained from the experimental data, 

is 0.375. In a coaxial swirling jet, the swirl number is defined as S = Jo P 1 ™ 1 ^ 

<pu 2 rdr ’ 

which U is the mean axial velocity W is the tangential mean velocity, and R is the annular 
jet radius. The inlet boundary conditions were set at 5mm downstream of the jet exit, at 
which the measured quantities in terms of mean and turbulence quantities are available. The 
results are obtained on a 81x81 non-uniform mesh with refinement in the recirculation 
regions and the entrance region. Grid independence was confirmed by performing another 
calculations with 101x101 grid cells. The difference between the two results is within 2%. 
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Figure 6 shows comparisons of the axial velocity along the centerline using the ASM 
closure and the two-equation model. In terms of the strength and location of the centerline 
reversal velocity, the ASM model yields better agreement than the k-e model. 
Comparisons of the predicted mean axial and tangential velocity profiles with experimental 
data are presented in Figure 7a and 7b. Both models show favorable agreements in the axial 
velocity profiles, but slight deviations in the radial components downstream. The good 
predictions obtained by the k-e model, which are in contrast to previous results in the 
literature, are mainly related to the second order difference scheme used in the present 
MAST code, and partly due to the relatively low swirl number situation. However, both 
models predict the rapid decay of tangential velocity to a solid-body rotation at the 
downstream locations. The cauculated tangential velocity recovers too quickly to the 
ultimate forced-vortex structure of confined swirling flows. This points to the deficiency of 
the ASM assumptions in which the “advection” of Reynolds stresses responds too quickly 
to the turbulent kinetic energy. The full differential RSM probably has to be used to 
correctly account for the history of turbulent momentum transport, which ultimately may 
slow down the recovery of the swirl velocity [35]. 

Comparisons of turbulent intensities and Reynolds stresses are shown in Figure 8a and 
8b. The predicted results of the two models qualitatively follow the trend of experimental 
data. In terms of the detailed profiles of turbulence properties, the ASM model conforms 
fairly well to the experimental data, while the k-e model predicts some relatively large 
deviations, especially, in the upstream region, where the non-isotropic turbulence prevails. 
The basic shortcomings of the k-e model for which the isotropic eddy viscosity 
assumption is employed indicates the inability to redistribute the Reynolds stresses. 

SSME Nozzle 

The geometry and chamber conditions follow the specifications used in [33]. The flow is 
subsonic at the inlet and supersonic at the nozzle oudet. A 81x71 grid with a very fine grid 
cluster near the nozzle wall was used to resolve detailed boundary layer flow structures. 
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Both the k - e/two-layer and the ASM/two-layer model were used for comparison purposes. 
Figure 9(a,b) show the iso-Mach lines of the non-reacting flow fields inside the nozzle. 
Figure 10(a,b) show the temperature contours. The difference between these two 
predictions is very small in terms of overall flow properties. The specific impulse of 100% 
power level is calculated as 513.53 seconds using the k -e/ model and is 513.79 using the 
ASM. The noted difference is in the predictions of turbulent quantities, such as the kinetic 
energy level, plotted in Figure 1 1 , for the near wall region at the nozzle exit. The different 
kinetic energy levels may have some impact on the near wall turbulent structure at the 
nozzle exit. 

This aspect can be investigated by carrying out a detailed calculation for the exit 
flow around the nozzle exit manifold. Geometry of the nozzle manifold is shown in Figure 
12 [34], Detailed temperature profiles along the nozzle wall was also supplied by Mr. 
Gross. The grids used for this calculation are shown in Figure 13(a) and (b). The 
calculated flow field characteristics from a previous whole SSME nozzle solution were 
used as inlet boundary conditions for the domain of interest. Specifically, the calculated 
profiles at the axial location x=117 [inches] from the nozzle throat were used. The 
boundary conditions away from wall were selected to ensure supersonic boundary 
conditions.. At the outer boundary for the calculation domain, either specified total pressure 
(for subsonic regions) or extrapolated boundary conditions (for supersonic regions) are 
specified. Two-layer treatment of the turbulence quantities as well as non-slip, isothermal 
boundary conditions were used at the wall. 

A 81x51 grid with about 15 grid points covering the near- wall inertia sublayer (y+ 
<50) was used for the calculations. The external total pressure was set to 1 Atm. The 
calculated flow fields, using the ASM model, were shown in Figure 14,15,16 in terms of 
Mach number contours, pressure contours, and vector plots. It is evident that a shock- 
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boundary layer interaction is experienced near the exit lip of the changing wall geometry. 
Due to the entrainment of external air, two recirculation bubbles are formed behind the 
oblique shock. A supersonic pocket region is also observed near the manifol exit due to the 
entrainment. Static pressures along the wall are plotted in Figure 17(a) and (b) using the 
two different turbulence models. Effects of wall temperatures are also shown. It can be 
seen that the k - e model moves the shock location slightly toward the nozzle exit and 
produces a lower pressure after the jump. The lower wall temperature also produced the 
same effect as the k - e by comparison to the ASM model. 

To investigate the effect of chamber pressure, another calculation was performed 
with a reduction in chamber pressure by 25% . From the pressure contour shown in 
Figure 1 8, the shock moved further inside the nozzle. The boundary layer around the shock 
region is still very thin. No analytical flow separation was observed. The final run 
involved molecular viscosity only, to simulate laminar flow. The boundary layer thickened 
and pushed the shock location upstream significantly as seen in Figure 19. However, still 
no flow separation was observed ahead of the shock formation. 

5. CONCLUSIONS AND RECOMMENDATIONS 

In this study, an Algebraic Stress Model (ASM) coupled with a two- layer near-wall 
treatment was developed and successfully implemented into the MAST code. To validate 
the model and the numerical implementation, a series of test case ranging from 
compressible flat plate flows, a recirculating flow, a confined swirling flow were computed 
and the results were compared with the base-line k - e model and available experimental 
data. Further applications including the entire SSME nozzle flow and the SSME exit wall 
flow around the nozzle manifold were studied. For recirculating flows and swirling flows, 
the ASM model shows improved results, compared to the k - e model, due to its ability to 
account for the non-isotropic effects. For SSME nozzle flow and exit flows around the 
nozzle manifold, the effects of second-order turbulence closure do not show significant 
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difference from the two-equation model calculations. The flow fields are dominated by 
shock/boundary layer interactions coupled with air entrainments from the outer boundaries. 

There are few suggestions for future work. First, the full Reynolds Stress Model 
(RSM) should be incorporated. One of the key difficulties in implementing the RSM is the 
specification of boundary conditions for Reynolds stresses and the near wall damping 
models. The two-layer model used in this study coupled with the ASM can be readily 
extended for RSM implementations. The RSM accounts for the history as well as transport 
of second-order stresses, which may better resolve detailed chock- boundary layer 
interactions. Second, the second-order closures should be implemented in three- 
dimensional numerical models to better resolve three dimensional effects. Third, the second 
order closure in terms of passive scalar transport such as turbulent heat fluxes (correlation 
of velocity fluctuation and velocity fluctuation) and turbulent mass fluxes (correlation 
between concentration fluctuation and velocity fluctuation) should be developed to better 
simulate mixing and combustion processes in thrust chamber. 
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Appendix 2 
Sample Inputs 


The MAST family computer programs consists of a set of subroutines controlled 
by a short main program. The fundamental structure can be found in the MAST user’s 
manual version 1.0 [37]. The updated capabilities, resulting from a previous study, were 
summarized in [38] for the version 1.1. Sample inputs for calculations of SSME thrust 
chamber flows and nozzle outlet manifold flows are given in Table A.l and A. 2. To 
activate the usage of the Algebraic Stress Model, keyword ASM is added in the 
TURBULEN Block as seen from Table A. 1 and A.2. 
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CONTROL COMPRES NCRT 3 OMGM 1 NCGM 50 

OMGD 1.0 PHI -1 OMGPHI 0.00 OMGT 0.50 OMGF 1.00 
ERRCG 1.0E-2 ERRM 1.0E-5 IMON 81 JMON 10 MONU 

; RESTART 

GRID NX 81 NY 71 AXISYM READXY 
BOUND 


1ST 

1 

I END 

1 

JST 

1 

JEND 

71 

INLET IPBC 3 

1ST 

1 

I END 

81 

JST 

71 

JEND 

71 

WALL U 0 V 0 TK 0 IPBC 3 

1ST 

1 

I END 

81 

JST 

1 

JEND 

1 

SYMMETRY IPBC 3 

1ST 

81 

I END 

81 

JST 

1 

JEND 

71 

OUTLET IPBC 3 

TURBULENT 


TKIN 3. 

TEIN 

1.E4 


ASM 



PROPERTY VISCOS -1 ; CALCULATE BY SUTHERLAND'S LAW 

PSTAG 20240946.90 TSTAG 3637. GAMMA 1.2 GMW 10.18 
SOLV U VP TEMP TK TE 

RUN DT l.E-7 DTMIN 2 . E-7 DTMAX 2 . E-5 CFLN 1.00 NSTEP 4 000 
NPR1 1 NPR2 100 NEX 18 

ENDJOB 


Table A.l. Input file of the whole SSME calculation using ASM model 


CONTROL COMPRES NCRT 3 OMGM 0 NCGM 50 

OMGD 1.0 PHI -1 OMGPHI 0.30 OMGT 0.30 OMGF 1.00 
ERRCG 1.0E-2 ERRM 3 . OE-5 IMON 81 JMON 10 MONU 

; RESTART 

GRID NX 81 NY 51 AXISYM READXY 
BOUND 

1ST 1 I END 1 JST 1 JEND 51 INLET IPBC 1 

1ST 1 I END 81 JST 51 JEND 51 WALL U 0 V 0 TK 0 IPBC 3 TEMP 300. 

1ST 1 I END 81 JST 1 JEND 1 OUTLET IPBC 2 

1ST 81 I END 81 JST 1 JEND 51 OUTLET IPBC 2 

TURBULENT TKIN 3. TEIN 1 . E4 ASM 

PROPERTY PREMAX 1.2E+5 PBACK l.E+5 

VISCOS -1 ; VISCOSITY IS CALCULATED BY SUTHERLAND'S LAW 
PIN l.E+5 TIN 300. GAMMA 1.2 GMW 10.18 
SOLV U VP TEMP TK TE 

RUN DT 5.00E-7 DTMIN 5. E-7 DTMAX 5. E-5 CFLN 1.00 NSTEP 4500 
NPR1 10 NPR2 100 NEX 20 

ENDJOB 


Table A. 2. Input file of the SSME exit flow using ASM model 
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7^/o with Mqo for adiabatic wall boundary condition. 
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m/ law 

Fig. 2 Variation of Cf/C/o with T W /T, IU , for M 
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plots of ?/ 4 t for adiabatic: wall boundary condition. 
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Fig. 4a Behaviors of k + and u'v' across the boundary 
layer for adiabatic wall boundary condition. 
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Fig. 4b Behaviors of *+ and uV across the boundary 
iayer for cold wall boundary condition ( T w /T uw = 0.4). 
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Fig. 5 Reynolds stress profiles for the backward-facing step turbulent 
flow (9:1), with data from [31] 
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ig. 7b Radial profiles of mean tangential velocity 
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Fig. 11 Behavior of near wall TKE at nozzle exit. 
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Fig. 13(a), Grid configurations for SSME nozzle exit manifold 
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Fig. 13(b), Close-up grids for Figure 13(a) 

43 













distance from sub— domain inlet along the wall, inches 


Fig. 17(a), Pressure levels along the wall near the nozzle exit using the 
ASM and k - e models 



distance from sub-domain inlet along the wall, inches 


Fig. 17(b), Effects of wall temperature on the wall pressure 
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